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ABSTRACT 


We  derive  explicit  formulas  for  the  second  moments  of  the 
absorption  time  matrices  in  the  Markov  Renewal  Branching  process 
These  formulas  may  easily  be  computationally  implemented  and  are 
useful  in  the  iterative  computation  of  the  semi-Markov  matrices, 
which  give  the  distributions  of  the  duration  and  of  the  number 
of  customers  served  in  a busy  period  of  a great  variety  of 
complex  queueing  models. 

KEY  WORDS 
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I.  Introduction 


This  note  is  a sequel  to  [1],  where  we  discussed  the  moments 
of  semi-Markov  matrices  associated  with  the  time  till  extinction 
in  a Markov  renewal  branching  process.  We  refer  to  [1]  for  the 
motivation  and  a general  description  of  the  problem  under 
discussion.  In  queueing  theory,  the  semi-Markov  matrices  of 
interest  are  related  to  the  number  of  customers  served  and  to 
the  time  duration  of  the  busy  period.  See  e.g.  [2]. 

Recalling  only  the  most  essential  definitions  from  [1],  we 
consider  a sequence  { A n ( - ) } of  substochastic  semi-Markov  matrices, 
with  Laplace-Stieltjes  transforms  (A* ( s ) } , such  that  the  sum 

oo 

(1)  A ( x ) = l A ( x ) , x >0 , 

n = 0 n 

is  an  irreducible,  regular  stochastic  semi-Markov  matrix.  The 
invariant  probability  vector  of  A = A(+«)  is  denoted  by  and  the 
vector  is  defined  by 

oo 

(2)  b = E nA  e , 

n=  1 n 

where  An  = An(+«),  for  n>0,  and  e = ( 1 , 1 , . . . , 1 ) 1 . 

Throughout  this  paper,  we  shall  assume  that  the  i rreduc i bi 1 i ty 
assumptions,  stated  in  [1],  hold  and  also  that 

• 

( 3 ) P = TT  6 < 1 . 

There  then  exists  a unique  sequence  of  substochastic  semi- 
Markov  matrices  (Gp(  •)},  such  that  the 

I 


2 


(4)  6* ( z , s ) = E z"  /"  e" sx  dG  (x), 

n = l 0 n 

defined  for  0<z<l  , s>0,  satisfies  the  nonlinear  matrix  equation 

(5)  G* ( z , s ) = zE  A* ( s ) [G* ( z , s ) ] . 

n = 0 n 

oo 

Moreover  the  matrix  G=G*(1,0)=  E G (+»)»  is  a positive 

n = l n 

stochastic  matrix,  which  satisfies 


(6)  G = E AG 
n = 0 n 


n 


The  invariant  probability  vector  of  G is  denoted  by  c[  and 
the  matrix  G is  defined  by  G..  = g.,  for  l<i,j<m.  The  matrix  I-G+G 

I J J 

is  known  to  be  nonsingular.  The  matrices  M and  N are  defined  by 


(7)  M = / xdE  G ( x ) = - [ TT  G* ( 1 , s ) ] , 

0 n=l  n ds  s=0 


N = E nG  (+«)  = [H_  G* ( z , 0 ) ] , 
n=l  n az  z=l 


and  pj=Me,  v/=Ne. 

For  future  reference,  we  define 


(8)  B ( 1 ) = E nA  , 
n = l n 

oo 

B ( 2 ) = E n (n-1 )A  , 
n = 2 n 


c ( 1 ) = / xdA(x)--[g—  A* ( s ) ] , 

0 05  s = 0 


£=B(1 )e, 

£2  = B ( 2 ) e_, 

3*  = C ( 1 )e. 


C(2)  = /VdA(x)  = [^-r  AMs)], 


ds 


s = 0 


££=C(2)e, 


D = E n/°°xdA  (x)  = - E n A*  ' ( 0+ ) , 6 = De, 


n = l 0 


n = l 


3 


E ( 1 ) =-  E A* 1 ( 0+ ) Gn  , 
n = 0 

00 

E ( 2 ) = Z A*"  ( 0+ ) Gn  . 
n = 0 n 


These  matrices  will  be  assumed  to  be  known  and  finite.  We  note 
that  if  G is  known  accurately,  then  E(l)  and  E(2)  can  in  principle 
be  computed,  with  E(l)e=8*  and  E(2)e=£^,  serving  as  accuracy  checks. 

By  routine  differentiations  in  Equation  (5)  and  using  results 
in  [1],  the  matrices  M and  N are  the  unique  solutions,  respectively 
to  the  matrix  equations 

(9)  M = E ( 1 ) + £ A nj;1GvMGn"1  "v, 
n= 1 n v=0 

00  n - 1 _ , 

N = G+r  A Z GVNG  'v, 
n= 1 n v=0 


and  the  vectors  and  v are  given  explicitly  by 


(10)  u = (I-G+G)[I-A+G-a(3)G]“V, 

v.  = ( I -G  + G ) [ I - A + G-  A ( fl)  G]  ^e  , 
where  a( p ) = di ag ( B1  , . . . , Bm) . 


Powerful  accuracy  checks  on  numerical  computations  are 
provided  by  the  formulas 

(11)  gu  = ( 1 - P ) “ 1 7T  b * , 2v  = (1-p)'1. 

The  purt  ^ of  this  paper  is  as  follows.  I*  '.ussing  the 

steady-state  probability  distributions  of  a variety  of  queueing 
models,  it  turns  out  that  only  the  vectors  £,  £ and  v and  the 
matrix  G are  of  essential  computational  Importance.  Even  for 
matrices  Ap  of  fairly  high  order, 


the  amount  of  comoutation  is 
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smal 1 to  moderate . 

In  contrast,  the  numerical  computation  of  the  semi-Markov 

co 

matrix  G(x)  = E G n ( x ) , for  x>0,  and  of  the  sequence  {Gn(+“)1  is 
n = 1 

a major  task.  The  most  promising  algorithms  involve  iterations 
of  the  nonlinear  operator,  corresponding  to  the  inverted  form  of 
Equation  (5). 

Such  algorithms  require  however  that  a priori  cut-off  points 
on  the  tails  of  the  probability  distributions  be  computed.  With 
information  on  the  second  moments  of  G(*)  and  {Gn(°°)>  available, 
we  can  limit  computation  to  terms  up  to  three  or  four  standard 
deviations  beyond  the  mean.  We  the  re fore  proceed  to  a discussion 
of  the  matrices  and  N2,  defined  by 

(12)  M2  = /°°x2dG(x)  = G*  ( 1 , s ) ] , 

0 ds  s=0 

00  - j 2 

N„  = E n ( n - 1 )G  (+«)  = [%.  G*(z,0)], 

6 n=2  n dz  z=l 

and  the  vectors  u_2  = M2i’  ancl  ^2  = N2-' 

Although  the  derivation  of  explicit  formulas  for  the  latter 
is  primarily  a matter  of  calculation,  the  intermediate  steps  are 
sufficiently  involved  that  it  appears  worthwhile  to  have  the 
resulting  expressions  available  in  the  literature. 

II.  The  Matrices  M^  and  N2. 

By  routine  differentiations  in  Formula  (5),  we  obtain  the 
equations 


5 


(13)  M9-  E A V GvM9Gn"1_v  = 

2 n -n  2 

n= 1 v=0 

°°  n-l  , 

E ( 2 ) - 2 E A*' (0+)  E GvMGn~  ~v 
n= 1 n v=0 

°°  n - 2 v 9 

+ 2 E A E E GrMG  MGrW-V, 

n=2  v=0  r=0 

and 

(14)  N9-  E A "e^  GvN9Gn‘1_v  = 

2 , n „ 2 

n=l  v=0 


2N-2G+2  E A E E Gr NGV‘ rNGn ~2 "v . 
n=2  nv=0  r=0 


For  practical  purposes,  the  explicit  expressions,  given  below, 
for  ^ anc*  ^.2  be  su  ffi  c i ent . The  numerical  solution  of 

Equations  (13)  and  (14)  is  feasible,  whenever  the  right-hand  sides 
can  be  computed  accurately.  The  following  results  convey  useful 
information  on  such  computation. 

Lemma  1 


For  any  square  matrix  X,  we  have 
o n - 2 v p 

(15)  lim  E E G XG  XG  = GXGXG 

n-voo  v = 0r  = 0 


(£Xe)2G. 


Proof 

n-l  n-l - 

Writing  U = E GvXGn“  _v,  and  noting  that  by  interchanging 
n v = 0 

the  order  of  the  summations,  we  obtain 

(16)  T * E E GrXGv'rXGn‘2‘v  = VgvXU  ,, 

n v=0  r=0  v=0  n_v"1 

so  that 

(17)  n-2T„  - n'1  "i*  6v*(n-v)-1Un.w_1  . 


In  [1],  Thm  4,  we  proved  that  n Un->-(<jXe)G.  Using  this 
result  and  repeating  verbatim  the  same  proof  for  the  expression 
in  (17),  the  stated  result  follows. 

Lemma  1 shows  that  the  truncation  error  at  the  index  n=K 

°°  n ~ 2 v r n ? 

in  the  infinite  sum  Z A z Z G MGV'  MGn"  v,  is  approximately 

n = 2 n v=0  r=0 

equal  to 

(18)  (l-p)'2Ue,*)2  E n2A  . 

n = K+l  n 

Similar  error  estimates  apply  to  the  other  infinite  series. 

The  matrices  Un  and  Tn  can  be  computed  recursively  for  use  in 
the  computation  of  the  right  hand  sides.  We  have  the  recurrence 
relations : 

(19)  U^X,  Un  = XGn"1+GUn  l , for  n>2, 

T2=X2,  Tn=XUn_2+GTn_i * for  n>3. 

III.  The  Vectors  p^  and  I2 ' 

Theorem  1 

The  vectors  p^  and  are  given  by 

(20)  E2  = 2[(3lL)I-M](I-G+G)"1E  + 

(I-G+G)[I-A+G-a(bJG]  ( 9 tL ) j®.+  (fly.) 2 0.2 

+2[C(1)+(9p)B(1 )-(3p)I](I-G+G)_1p}, 
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(21)  v2  = 2[(1-p)_1 I-N](I-G+G)"1v  + 2(l-p)"1pv  + 

(I-G  + G)[I-A+G-A(6)G]_1 { (1  -p)'2£2  + 

2(1-p)'1[B(1  J-pIKI-G+G)"1^}. 

Proof 

The  proofs  of  both  formulas  are  similar,  so  we  shall  only 
give  the  main  computational  steps  of  the  proof  of  (20).  Several 
non-obvious  steps  are  involved. 

Multiplying  the  left  hand  side  of  (13)  on  the  right  by  e, 
we  obtain: 

<»  n - 1 __  _ , 

(22)  s A £ GVu_2  = + )G](I-G+G)  u_2 

n=l  v=0 

= [I-A+G-a(6)G] ( r -G+G ) ~ 1 u2* 

Multiplying  the  terms  in  the  right  hand  side  by  e,  the  first  term 
clearly  yields  . The  second  term  yields 
® n - 1 

(23)  2 I A*’ (0+)  E G\  = 

n = l n v = 0 

00  oo 

2[A* 1 ( 0+  ) - E A* 1 ( 0+ ) Gn+  E nA* 1 ( 0+  ) G] ( I -G  + G ) ~ y 
n=0  n n=l  n 

=2[E(1 )-C(l )](I-G+G)"]u-2(2u)6. 


The  third  term  is  more  involved.  We  obtain 

(24)  2 E A z E GrMGV'rp  = 

n= 2 n v=0  rr0 

2 E A % I GrM(GV'r-GV+1"r  + G)(I-G  + G)‘V  = 

n=2  v=0  r=0 

00  n-1  , ® n- 1 , 

2 E An  E GVM ( I -G  + G ) " y -2  E A E G vMGn ' 1 " v ( I - G+G ) 
n=l  n v=0  n= 1 v*0 

o»  n - 2 v 
+2  E A e £ G MGu. 
n=2  v=0  r=0 
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The  first  two  of  the  latter  terms  simplify  to 

(25)  2[A-G+A(g)G] ( I-G+G)-1M( I-G+G)_1u 

-2[M-E(1 )](I-G+G)_1p  = 2E(1 ) ( I -G+G) _1 y 
-2[I-A+G-A(e)G](I-G+G)"'IM(I-G+6)“1u . 


The  last  term  in  (24)  yields 
00  n-2  v 

(26)  2(gu)  T A 1 1 6 w = 

n=2  n v=0  r=0 


2(fly.)  I A E [I-GV  + 1 + (v+l  )G]  ( I-G  + G)-1p  , 
n = 2 n v = 0 


a.  n-2  i 00  » n-1 

(27)  E A E [ I -Gv  +( v+1  )G]  = I nA  ■ E A E Gv 

n = 2 n v = 0 n = 2 n n = 2 n v = 0 

+i  E n (n-1  )A  G = B ( 1 ) - [A-G  + a ( 0 ) G]  ( I -G  + G ) _1 +L  ( 0?  ) G , 
c n=2  n c c 

so  that  the  expression  in  (26)  becomes 

(28)  (9y.)2B2+2(2p)[B(l  )-I](I-G+G)_1u 
+2(gp)[I-A+G-A(0)G](I-G+G)‘2p. 

After  collecting  terms  and  multiplying  by  the  inverse  of  the 
coefficient  matrix  in  Formula  (22),  we  obtain  the  explicit 
formula  (20). 

Corollary  1 

In  the  M/SM/1  queue  with  arrival  rate  X and  service  time  semi- 
Markov  matrix  A(*)»  we  have  the  simplified  formulas 
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(29)  p = XTre*,  B(1)=AC(1),  b(2)  = a2c(2), 

£=A|3*,  B2=A2£*  i=xi2 

gp  = A 1 p ( 1 -p  ) 1 , 

w.2  = 2[A"1p(l-p)"1I-M](I-G  + G)“1iL  + 

( I-G  + G)[I-A  + G-AA(g*)G]_1  { ( 1 - p ) “ 2 6_|  + 
2(l-p)"1[C(l  ) - A ' 1 p I ] ( I -G  + G ) " 1 vl)  , 

y_2  = 2[(l-p)  ^ I— N](T— G+G)  ^v+2(l-p)  1 pv  + 

( I -G+G) [I-A+G-AA(b*)G]_1 { (1  - p)'2A2£*  + 
2(l-p)_1[XC(l )-PI](I-G+G)-1v}. 


Remark 

In  the  scalar  case  (M/G/l),  where  A = G = G = 1 , AB*=p,  v=(l-p)  \ 
p = A ^p(l-p)-^,  the  latter  two  formulas  reduce  to  the  classical 
f ormu las. 

U2=(l-p)  ^&2>  v^O-p)  3 A2 B*  + 2p ( 1 -p ) 2. 

IV.  Numerical  Computation 

Assuming  the  required  moment  matrices,  other  than  E(l)  and 
E(2),  defined  in  (8)  are  known,  the  first  step  is  the  numerical 
computation  of  G.  As  discussed  in  [1],  we  found  that  successive 
substitutions  in  the  equation 

(30)  G = ( I -A , )-1  5:  A Gv, 

1 v = 0 v 

v/1 

starting  with  G=0,  converge  rapidly  except  for  p very  close  to  one. 

The  evaluation  of  G and  E(l)  is  routine.  Iv  is  not  necessary 
to  compute  the  inverse  of  (I-G+G),  since  only  the  vectors 
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(I-G+G)  \ and  (I-G  + G)”^v.  are  needed  and  by  Formula  (10),  it  is 
clear  that  these  can  be  evaluated  without  inverting  the  matrix 
I-G+G . 

The  main  computational  effort  goes  into  the  evaluation  of 

the  matrices  M and  N.  As  shown  in  [1],  the  equations  (9)  are  each 

2 2 

equivalent  to  systems  of  m linear  equations  in  m unknowns. 

Except  for  very  small  values  of  m,  it  is  impractical  however  to 
compute  the  required  coefficient  matrix. 

We  have  successfully  computed  matrices  M and  N for  m as  large 
as  thirty,  without  expending  excessive  computer  times  by  straight- 
forward successive  substitutions  and  by  using  the  first  of  the 
recurrence  relations  in  (19)  to  compute  the  required  matrices  Un . 

The  speed  of  convergence  can  be  substantially  improved  by 
writing  the  equations  (9)  in  the  form 

(31)  X=  [ I - Z A Gn" 1 ] 1[Y+  Z A "zVxg"'1^], 
n-l  n = 2 v = u 

with  X=M,  Y=E(1)  and  X=N,  Y=G,  respectively. 

Ac  know! edqemen  t 

The  author  thanks  Mr.  D.  Lucantoni  for  checking  the  tedious 
calculations  involved  in  this  paper. 
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